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Abstract 



At high densities fluids of strongly dipolar spherical particles exhibit spon- 
taneous long-ranged orientational order. Typically, due to demagnetization ef- 
fects induced by the long range of the dipolar interactions, the magnetization 
structure is spatially inhomogeneous and depends on the shape of the sample. 
We determine this structure for a cubic sample by the free minimization of 
an appropriate microscopic density functional using simulated annealing. We 
find a vortex structure resembling four domains separated by four domain 
walls whose thickness increases proportional to the system size L. There are 
indications that for large L the whole configuration scales with the system 
size. Near the axis of the mainly planar vortex structure the direction of the 
magnetization escapes into the third dimension or, at higher temperatures, 
the absolute value of the magnetization is strongly reduced. Thus the orienta- 
tional order is characterized by two point defects at the top and the bottom 
of the sample, respectively. The equilibrium structure in an external field and 
the transition to a homogeneous magnetization for strong fields are analyzed, 
too. 

PACS numbers: 75.50.Mm, 75.60.Kw, 61.30.Cz 
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I. INTRODUCTION 



At high densities strongly dipolar fluids can exhibit an orient at ionally ordered liquid 
phase characterized by a spontaneous magnetization. This has been demonstrated for the 
first time by Monte Carlo simulations of dipolar soft spheres by Wei and Patey and has 
been confirmed by simulations of hard spheres and of Stockmayer fluids . Theoretical 
approaches leading to the same conclusion include a mean-field theory for a dipolar lattice 
gas 1]^, a generalized van der Waals theory and different kinds of density- functional 
theory [p|-pT[]. At low densities the typical configurations of dipolar fluids exhibit chain 
formation |T^-|14|,^ which may inhibit the phase transition from the isotropic vapor phase 
to the isotropic liquid phase [p!5|-[T8[|. If the dipole moments are electric as in molecular 
liquids the orientationally ordered phase exhibits a spontaneous polarization; in the case of 
magnetic dipoles as in ferrofluids, i.e., colloidal suspensions of permanently ferromagnetic 
particles, one speaks of ferromagnetic order and a spontaneous magnetization. In this paper 
we adopt the magnetic language, keeping in mind that completely analogous phenomena 
occur in the electric case (as long as no free charge carriers are present in the fluid). 

In both simulations and analytic theories the dipolar forces must be treated carefully 
due to their long range which may give rise to effects depending on the shape of the sample. 
It turns out that for all sample shapes with the exception of a long needle the equilib- 
rium configuration is inhomogeneous with a spatially varying magnetization M(r) |lT9| , pO| . 
This leads to a shape independent free energy, as is expected on general grounds The 
highly non-trivial problem to determine explicitly the spatial distribution of this inhomo- 
geneous magnetization for a given sample shape has not yet been solved satisfactorily. In 
simulations usually a homogeneous magnetization is enforced by using an infinitely per- 
meable surrounding of the periodically repeated simulation cells. If instead the sample is 
surrounded by vacuum the simulation cell splits into two domains with opposing magnetiza- 
tions 0. But this structure is clearly induced by the artificial periodic boundary conditions. 
We addressed this problem recently |^ for the experimentally more relevant case of open 
boundaries and obtained the following general characterization of the equilibrium configura- 
tions: On a macroscopic scale the absolute value of M(r) is constant, its divergence is zero, 
and at the surfaces the normal component vanishes. However, these properties do not yet 
enable one to construct the configuration for a given shape. For a cubic sample we deter- 
mined [0 the most stable configuration under the constraint of sharp domain boundaries. 
It consists of four triangular domains with 90° domain walls in between. In contrast to solid 
ferromagnets the number of domains does not increase with the system size, but remains 
as small as possible. On the other hand the assumption of sharp domain boundaries is not 
justified for liquid systems. Due to the lack of a lattice anisotropy one rather expects very 
thick domain walls p2 |. 

Therefore, in order to analyze the domain configuration in more detail in the present 
work we have performed a numerical minimization of an approximate microscopic density 
functional which has been used before [|I^,|r2,^] to describe the ferromagnetically ordered 
fluid. We have made no a priori assumptions regarding the symmetry of the domain structure 
and have minimized with respect to a large number (10^-10^) of parameters which repre- 
sent the local magnetization at mesh points within the sample. This approach is similar in 
spirit to the determination of the magnetization structure of solid ferromagnetic particles of 
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micrometer size in the framework of micromagnetic theory [p^-|25[]. But there are two main 
differences. First, we do not constrain the absolute value of the magnetization, thus allowing 
for the formation of less ordered regions. Second, we work on the microscopic scale deter- 
mined by the particle diameter, which enables us to examine in detail the overall structure 
as well as the structure of domain walls and the cores of topological defects. 

A brief account of some of our results has been published in Ref. p6|]. Here we provide the 
important derivation of the functional from the underlying model, discuss the implementa- 
tion of the simulated annealing, and give a thorough analysis of the resulting configurations. 
We also analyze the consequences of the necessary discretization and of the finite sample 
size. Furthermore we discuss in detail the orientational structure in an external field, which 
is relevant for the experiments with the recently discovered metallic liquid ferromagnets 

m. 



II. MODEL AND THEORETICAL APPROACH 

A. Density-functional theory for Stockmayer fluids 

We consider Stockmayer fluids consisting of spherical particles with fixed embedded point 
dipoles which interact via pairwise dispersion and dipolar forces. The interaction potential 
10 = 10^ + Wdip is the sum of the Lennard- Jones potential 

and the dipolar potential 

2 

771 

Wdip{ri2,uj,uj') = -Y- (m(t^) ■ m(u;') - 3{m{uj) ■ ri2)(m(u;') ■ ri2)) 6(ri2 - a). (2) 



12 



ri2 = r — r' = ri2ri2 denotes the interparticle vector, uj = {6, (p) and uj' are the orientations 
of the dipole moments at r and r', respectively, with an absolute value m; hats indicate unit 
vectors. At short distances the interaction is cut off by the Heaviside function B. 

In order to study the spatially inhomogeneous configurations within a ferromagnetically 



ordered fluid we employ the density-functional theory introduced in Ref. and for which 
we have worked out analytical results previously |]TD|,|TT|,^ . 

The configurations of the fluid are described by the spatially constant number density 
p and the normalized orientational distribution a{r,uj) so that the probability density for 
finding a particle at point r with the orientation uj is p{r,uj) = pa{r,u). The free energy 
density functional is given by 

F{p, [a{r, uj)]; T) = Vfnsip, T) + For + F,^ + Fh. (3) 

V is the sample volume and fus is the free energy density of the hard sphere reference 
system characterized by an effective temperature dependent radius PU|. The second term 
(/3 = ll{kBT)) 
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or 



J d^r J (icua(r,tu) ln[47ra(r, tu)] (4) 



takes into account the loss of entropy if the orientational distribution is not isotropic, i.e., 
different from l/(47r). Here and in the following the integrations over the angles u are taken 
over the unit sphere. Within the low-density approximation the excess contribution due to 



the long-ranged part of the interaction potential has the form |^ 

^''^ ^~2pjy jy / ^)"(^'' ^')0(n2 - ^)/(ri2, uo, uo') (5) 

with the Mayer function 

/(ri2,w,cj') = -1 + exp[-/5u7(ri2,u;,u;')]- (6) 

This function can be expanded in terms of the rotational invariants {C{lil2l,mim2m) are 
Clebsch-Gordan coefficients) 

^l^l^l{uJ,Uj',UJi2) = 5Z C{lil2l,mim2m)Yi^rn.A^^)yi2rn2i^')yinX^12) (7) 
mi,m2,m 

with coefficients / depending only on the distance between the particles: 

/(ri2, a;, u') = ^ fhi2iiri2)^hi2ii^^, uj\ uJi2)- (8) 

Finally, the interaction energy with a homogeneous external field H is 

Fh = —pmH J d^r J dcu a {r,uj) cos 'j (9) 

where 7 is the angle between u and the direction of the field. We take into account the 
external potential due to the container walls only summarily in that they provide the con- 
finement of the fiuid to V because we are interested in the bulk behavior of the fiuid. In this 
spirit we also do not consider spatial variations of the number density p(r) = J du p{r, uj) 
in the vicinity of the walls. Moreover, in Eqs. (0)-® we have assumed that p(r) is constant 
throughout the sample even if a{r,uj) varies deep inside the sample. We expect that due to 
the small compressibility of the fiuid in the orientationally ordered phase possible variations 
of p(r), e.g. inside the domain walls, are only small and thus not significant. 

In the next step the orientational distribution is expanded into spherical harmonics: 

00 I 

a{r,u) = ^ ^ fiim{r)Yim{uj) (10) 

1=0 m=-l 

with Poo = l/"\/47r due to the normahzation J duj a{r,uj) = 1 and pf^ = {—l)"^ni-m because 
a is real. Inserting Eqs. and (0) into Eq. (|^) one finds for the excess free energy 



4 



17111712171) 



h,h,l mi,m2,m 
'3^ / J3^', 
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(11) 



In order to make the free minimization of the density functional with respect to the set 
{fiim{T^)} numerically feasible we must refrain from determining the full orientational distri- 
bution a{r, u) at each point r. Instead we focus on the reduced information provided by the 
dimensionless local magnetization 



M(r) 



duj a{r, u)m{Lj) 
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(12) 



(The actual magnetization is A^(r) = ^/AttJs pmM{r).) This implies that in Eq. (|ToD only 
terms up to / = 1 have to be taken into account. Within this approximation those contribu- 
tions to the excess free energy which depend on the orientational order are given by 



2/3^12^ 
2/3 



d'r / ciVM(r)-M(r')e(ri2-a)/no(ri2) 
V Jv 



(13) 
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d'r / rfVM(r)t(ri2)M(r')e(ri2-fr)/ii2(ri2) 
V Jv 



with the tensor Tij(r 



Sij — ^fifj. We note that in contrast to Ref. |]20[ we do not have to 



separate long- and short-ranged contributions to fiu because within our numerical approach 
we are always dealing with finite systems without carrying out the thermodynamic limit. 



Analytic expressions for the functions fm are provided by Eqs. (B33) and (B34) in Ref. ||30 
from which one obtains the following expansions in terms of m^: 



/iio(r) 
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I __p-PwLj{r) 
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+ 0(m 



ION 



(14) 



f 



of 



To lowest order in m? the function fiuif) contains the characteristic dependence 
the dipolar potential whereas the function /no arises only due to the cubic and higher order 
terms in the expansion of the Mayer function. (For the actual calculations we have used the 
full expressions as given in Ref. |^| instead of the expansions in Eq. (p!^.) 

In the following we specificly analyze a cubic volume of side length L containing 
the fluid. The spatial integrations are discretized by introducing a (simple cubic) lat- 
tice with lattice constant a = L/N and lattice vectors R = a(ni, ria, 7^3), with ni G 
{-{N - l)/2, . . . , -1, 0, 1, . . . , (A^ - l)/2} {N is even), so that 



AF,, 



R R' 



WijiK 



R')Mj-(R'). 



(15) 
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The interaction tensor Wij is given by 



2A7T 



Q{R-a) 



I, J = 1,2,3, (16) 



so that R) = Wij{H). (The value of a /a is taken to be noninteger; otherwise there 

are distances between lattice points which correspond to the discontinuity of the Heaviside 
function.) An alternative formulation, which will prove to be helpful later on, is 



A^e. = |^a3$^5^if,(R)M,(R) 



R i 



with the local fields 



R' j 



(17) 



(18) 



The entropic term given by Eq. (^) can be simplified using the fact that by applying a 
suitable rotation the orientational distribution at a given point r can be cast into the form 



a{r,u;) = ^ + M(r) W ^ cos^ 



(19) 



with M(r) = |M(r)|. Due to the angular integration in Eq. (^ this rotation does not alter 
the value of Fnr so that 



For = ^ j d^r j duj{^ + M{r)^l^cose^ In (l + M(r)v^cos^) . 



(20) 



The expansion of the logarithm yields in the discrete version 

(v^M(R))2" 



R n= 



^ {2n - l)2n{2n + 1) 



(21) 



5^.(M(R)). 



R 



Finally, the term due to the external field, which is taken to point into the z direction, is 

pmHa"' 

R 

Thus the free energy difference F^f between the /erromagnetic and the zsotropic phases is 




(22) 



"i/ — For + ^Fex + Fff 



(23) 



where the individual terms are given by Eqs. ([T5|), (|2T|), and (|22D, respectively. 
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B. Simulated Annealing 



We have minimized the free energy in Eq. (^) with respect to the magnetization con- 
figuration {M(R)} by using the simulated anneahng algorithm In the first step the 



interaction tensor Wij is determined once for all relevant distances between the lattice sites 
within V = L^. In the second step the local fields -ffj(R) are calculated for an initial choice 
of the configuration M^*™)(R) according to Eq. (|18]). In each following step a lattice site 
Ro is chosen and a new magnetization M'(Ro) is proposed by changing each component 
i = 1, 2, 3 by a random value AMi = M;(Ro) - Mi(Ro) between -kM(Ro) and +/tM(Ro) 
= 0.1 turned out to be a suitable choice). The resulting change in free energy is 



s(M'(R„)) - s(M(R,)) + J2 H,{Tt<,)AM. 




Ayr 

—pmHa^AM^. (24) 



There is no term quadratic in AMi due to Wij(R = 0) = 0. The proposed change of 
the configuration is accepted with certainty if AFif is negative and with a probability 
exp{—AFif / {ksTs)) if it is positive. Tg is the control temperature of the annealing algo- 
rithm which is decreased slowly during the minimization. In case of acceptance the new field 
i?-(R) at each site is calculated according to (see Eq. (|ll 



H'i{K) = Hi{K) + a^Yl " Ro) AMj. (25) 

j 

The time required for this computational step, which is consuming most of the CPU time, 
is of the order 0{N'^). The advantage of using the local fields is that they need not be 
updated if the proposed change is rejected, which happens quite often especially at low 
control temperatures near the end of a run, while the acceptance decision itself can be 
reached very fast. 

The control temperature is lowered by a factor of 0.95 after 3A^^ successful changes or 
15A^^ trials. The algorithm is terminated if no successful step occurred during the 15 A^^ 
trials. With the assumption that the number of temperature steps is independent of N the 
total computing time should scale as N^. Actually we found a scaling exponent between 
6 and 7. Typical CPU times for one run on a DEC alpha workstation were 3.2 hours for 
N = 1Q and 25 hours for iV = 22. 

The minimizing configurations can be found starting from completely random initial 
states, but this requires a relatively large initial value of Tg and therefore very long runs. 
For this reason we started in almost all cases from a configuration which had been obtained 
as a minimum for other parameter values. It is partly randomized during the first phase of 
the algorithm by applying an appropriate initial control temperature so that a large fraction 
of proposals is accepted. By testing in some cases different starting configurations we took 
care not to bias the final result by a prejudiced initial guess. 



III. RESULTS AND DISCUSSION 
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A. Magnetization structure 



As our standard values of the thermodynamic parameters we chose m* = ^Jw? ja'^e = 
1.5, T* = ksT/e = 2.25, and p* = pa^ = 0.94. This leads to a thermodynamic state which, 
within the present density-functional theory approximation, lies deep in the ferromagnetic 
liquid phase. In view of the numerical challenges described in Subcsec. [11 B| we have examined 
system sizes L* = L/a between 4.8 and 12 using lattices consisting of = 10, 12, ... ,24 
sites in each dimension. (Although we used an ordinary workstation our maximum system 
size is even 30% larger than the system consisting of 22'^ sites examined by Williams and 
Dunlop [0] on a supercomputer in 1989.) 

The result of a minimization run is a three-component vector field M(R) representing the 
magnetization structure within the cubic volume. In order to visualize this field we display 
sections parallel to the faces of the cube with the magnetization projected onto the section 
plane. Since one component is lost due to this projection the absolute value |M(R)| cannot 
be inferred from these figures. We adopt a reference frame which has its origin in the center of 
the cube. Figure 0(a) shows a section perpendicular to the z axis a.t z* = z/a = 0.18 ~ 0, i.e. 
close to the center of the cube, for L* = 7.2 and = 20. The overall picture is that of a vortex 
of closed magnetization lines circulating around the z axis. In this context it is interesting to 
note that clusters of some ten to hundred dipolar particles also exhibit a vortex structure at 
low temperatures p2| , p3| . This structure leads to divH = for the resulting magnetic field 



H. A closer look at the structure in Fig. |I|(a) reveals that it may be described as composed 
of four domains with an approximately constant magnetization separated by broad domain 
walls along the diagonals of the square within which the direction of M changes continuously. 
This configuration resembles the triangular structure which has been found to be the most 
stable structure under the constraint of sharp domain walls [^. A similar structure has also 
been found to be the most favorable one in cubic magnetite particles just above the critical 
single-domain size, while more complicated structures occur in this case for larger particles 
2^,|2S|. In a section near the bottom of the cube {z* = —3.06 ~ —L*/2, Fig. 0(b)) the domain 



walls are slightly shifted off the diagonals into the direction opposite to the circulation. The 
reverse situation is found near the top. However, in contrast to the triangular structure 
studied in Ref. [^], the magnetization is not confined to a plane. As can be seen in the 



sections perpendicular to the y axis in Fig. |^ there is a nonvanishing z component, leading to 
an "escape into the third dimension" , which is particularly pronounced near the vortex axis 
(Figs. H(b) and (c)). This mechanism avoids the formation of a topological defect along the 
z axis. This is in accordance with general considerations showing that line singularities are 



topologically unstable in a system of three-component spins in three spatial dimensions [34 



which means that they can always be removed by continuous local modifications. However, 
near the top and the bottom face of the cube the z component decreases in order to avoid a 
large normal component at the surface which would produce an unfavorable demagnetization 



field |20]. Thus two topologically stable point singularities near the centers of the bottom 



and top surfaces are inevitable. A closer look at the magnetization structure reveals a strong 
decrease of the absolute value of M in the core of these defects (see Fig. H). The fact that 
the size of this less ordered core also scales with the system size underlines the importance of 
including the absolute value of the magnetization as a minimization parameter. If instead the 
assumption of constant magnitude of M were applied, as is usually done in the micromagnetic 
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theory, the defect could not be described correctly even on a mesoscopic scale. 

Very similar structures have been found for all values of L and A^. In Fig. ^ we compare 
the "degree of escape" for different system sizes with fixed = 20. The cosine of the polar 
angle of M (cos^ = M3/M) is plotted as a function of the distance r = from the 

center along the diagonals (i.e., the lines x = ±y with z = const) and the center lines (i.e., 
the lines x = and y = with z = const). With the exception of the smallest system the 
scahng property 

M(r/a, L/a) = M^") (r/L) (26) 

is approximately fulfilled, so that one can surmise that this holds also in the thermodynamic 
limit. (Since M is dimensionless, for a given thermodynamic state it can only depend on 
r/a and L/a.) The function M*^"^ represents the global texture in the thermodynamic limit. 
Near the edges of the cube, i.e., near the corners in the projection plane, there is an escape 
into the opposite direction, but combined with a strong decrease of the absolute value of 
the magnetization (see Fig. |^). Upon moving outwards from the center cos^ decreases faster 
along the diagonals than along the center lines, indicating that there is not a circular but 
rather a square symmetry. 

In the inner part of the sample the absolute value M(R) = |M(R)| (Fig. |^) is approx- 
imately constant and independent of L. It decreases near the surface and near the edges. 
This less ordered surface layer thickens relative to the system size as L decreases. Here, 
too, scaling with the system size for large L (Eq. (p6D ) is compatible with the data. In the 
smallest system we studied [L* = 4.8) M decreases also near the center which indicates a 
second mode for avoiding a line singularity besides the "escape into the third dimension". 
This mode appears also for the larger systems near the ferromagnetic-isotropic phase tran- 
sition and might be induced by the proximity of this transition for L* = 4.8 (see Sec. [Ill C| ). 
In Fig. ^ we present the average magnetization within each plane parallel to the xy plane 
{M)xy = N'"^ M{nia,n2a, z) as a function of the height z. Again there is a clear 

decrease near the surfaces while in the central region {M)xy is constant. The average {M)xy 
attains its thermodynamic limit L — 00 more rapidly than cos^ or M(R) (see Figs. ^ and 

For small L the surface disordered region shrinks on the scale of L but in the thermo- 
dynamic limit it remains proportional to L. For large L we find for the excess quantity 
j% dz [(M).,(0) - {M),y{z)] /(M).,(0) ^ 0.08L. 

As shown in Fig. ^ the minimum value of the free energy density f*^ = {Fif / L^){a^ / e) 
exhibits a relatively strong, oscillatory dependence on the lattice constant a = L/N with a 
decreasing amplitude for increasing values of A^. This figure reveals that minima occur when 
a/a = L*/N is close to the inverse of an integer, which indicates a strong discretization 
effect. For the larger values of L we could not reach the region where the oscillations have 
died out, which inhibits a further finite-size analysis. Similar oscillations were found in the 
spatially averaged value of M(R) and for L* = 4.8 also in the degree of the rotation out of 
the plane. The precision to which the minimum value of f*j can be determined by simulated 
annealing for a given set of parameters is much higher than the discretization effect described 
above; from the dependence f*j{Ts) we estimate the minimization error to be of the order 
of 0.01. 
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B. Domain walls 



In solid ferromagnets the walls between adjacent domains have a microscopic thickness 
(i.e., it does not scale with the system size) which is determined by the competition between 
the exchange energy resulting from the spin coupling and the anisotropy energy due to the 
lattice structure which causes easy axes for the magnetization. Since there is no such lattice 
anisotropy in liquid ferromagnets de Gennes and Pincus p2| surmised that consequently 
there are also no domain walls. Below we shall argue that in the case of cubic samples as 
considered here the thickness of the domain walls is proportional to the system size and 
thus diverges in the thermodynamic limit. Thus one is left with a question of terminology 
whether one still speaks of walls, but certainly the behavior is qualitatively different from 
that in solids. 

In order to analyze the properties of the fluid domain walls in the finite cube we con- 
sider the behavior of the dimensionless magnetization M along straight paths normal to the 
wall. Along these paths M changes continuously between the magnetization directions of the 
adjacent domains. In the present case these domains have the quasi-triangular structure indi- 
cated in Fig. |l](a) so that the orientational order between neighboring domains differs by an 
angle of 7r/2. Except near the vortex axis we find a Neel type of wall, i.e., the magnetization 
vector rotates mainly within the plane spanned by the asymptotic orientations deep inside 
the adjacent domains; in Fig. |l|(a) this is the xy plane. In contrast, in bulk solid ferromagnets 
one usually observes Bloch walls with the magnetization vector rotating out of plane on a 



cone around the wall normal only in thin solid films ||3^ and small particles Neel 
walls do occur, too. In Fig. ^ the relevant transversal component Mt^^^y = {Mi + M2)/\/2 is 
shown as a function of the normal coordinate n (compare Fig. |l|(a)). 

We define the wall thickness 5 as the distance between the points where the tangent 
to the curve at its inflection point reaches the extreme values of Mt^xy (see Fig. |^). The 
dependence of 5 on the distance r of the wall normal from the center of the cube is displayed 
in Fig. 1^ (all paths lie in the plane z = —a/2, i.e., close to the center of the cube). The 
wall thickness decreases near the edge of the cube where the normal paths hardly reach the 
region of homogeneous magnetization (see Fig. |I|). The thickness increases near the center 
r = due to the "escape" region (see Fig. 0(c)). However, it is approximately constant in 
the range r/L = 0.2. ..0.35 and there is no monotonic trend as function of the number of 
lattice sites. Its dependence on the system size L is analyzed in Fig. |10[ Here we selected 
the values of 5 at (or, due to the discrete lattice, close to) r/L = 1/2'^/^ corresponding to 
half the distance between the center and the edge. For each size L we display the results 
obtained for different values of which render an estimate of the uncertainty caused by the 
finite lattice constant a = L/N. A slight decrease of 6/L with increasing L indicates that 
the domains are getting sharper. However, most probably the data can be extrapolated to 
a finite limit of 6/L for L — > oo which would be in accordance with the proposed scaling 
behavior in Eq. (|26|). 



C. Temperature dependence and critical point 

Starting from the standard value T* = ksT/t = 2.25 used up to here we have increased 
the temperature at the fixed density p* = 0.94 in order to examine the structural changes 
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upon approaching the ferromagnetic-isotropic transition. We have chosen = 16 and stud- 
ied the system sizes L* = 4.8 and L* = 9.6. The spatially averaged absolute value of the 
magnetization decreases and finally vanishes at a temperature Tc{L) in accordance with 
the inherent mean field approximation (see Fig. |ll]). (We define the finite-size critical tem- 
perature as the limiting temperature above which no configurations with a negative free 
energy difference Fij are found by the minimization algorithm.) As expected this finite-size 
(strictly speaking quasi-) critical temperature is lower for smaller systems. From Eq. (7.10) 
in Ref. [11] we infer T*{L — > oo) = 3.04 for the parameters used here. The evolution of the 
magnetization structure is analyzed in Figs. O and O for L* = 9.6. The escape into the z 
direction near the vortex axis is strongly reduced at higher temperatures (Fig. |12D while the 
absolute value of the magnetization in this region decreases more rapidly than for interme- 
diate values of r/L (Fig. |T3|). Thus a column of less ordered fluid develops around the vortex 
axis. These effects are even more pronounced in the smaller system, as has been already 
suggested by Figs. ^ and ||. Furthermore Fig. ^ demonstrates that the surface layer with 
reduced orientational order thickens, which is a consequence of the increasing correlation 
length upon approaching the phase transition. The domain wall thickness as deflned in the 
preceding subsection also increases slightly but this change is smaller than the uncertainty 
caused by the flnite lattice constant. 



IV. EXTERNAL FIELD 

Up to now all results refer to zero external fleld. If an external fleld is applied the dipolar 
particles tend to align along the fleld direction. The resulting transition from the inhomo- 
geneous zero-fleld conflguration to the homogeneously magnetized state in the presence of 
strong flelds can also be examined within the present approach. 

The fleld destroys the equivalence of the three perpendicular directions of the cubic axes. 
(One should keep in mind that below Tc{L) this equivalence is also spontaneously broken in 
zero fleld.) We have applied the fleld normal to the surfaces of the cube either parallel or 
perpendicular to the spontaneously chosen vortex axis of the zero-fleld conflguration, which 
was used as an initial guess for the minimization algorithm. The most stable conflgurations 
were found when the external fleld is parallel to the vortex axis. (In principle, equivalent 
but rotated conflgurations should be obtained if at the beginning of the algorithm the fleld 
is applied perpendicular to this axis. Since, however, only a medium value for the initial 
control temperature was employed, for this latter choice of the initial guess the system could 
not reach the equilibrium structure which in this case differs signiflcantly from the initial 
conflguration.) The relative stability of the resulting conflgurations can be judged on the 



basis of the corresponding value of the free energy. A typical result is depicted in Figs. O 



and 15. The section parallel to the xy plane, i.e., perpendicular to the external fleld, exhibits 



smaller absolute values of the magnetization components orthogonal to the fleld direction 



than for H = (Fig. |T^), but the structure is very similar. On the other hand there is a 



substantial increase of the magnetization component parallel to the fleld as shown in the 



sections y = const in Fig. |T^. Thus the overall structure is essentially preserved, but at 
each point the magnetization is rotated into the fleld direction by a certain amount. This 
demonstrates that this behavior, which was anticipated in Ref. [^, corresponds indeed 
to the equilibrium conflguration. If the fleld is sufficiently strong the vortex structure is 
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lost at a critical field strength He- This is shown in Fig. ^ by the field dependence of 
the averaged parallel component (My) = XIr ^sl^-) of the angular component 
{M^) = iV"^ERe^(R) ■ M(R), where e^(R) = {R2, -Ri,0)/\Ti\. Below the parallel 
component increases approximately linearly while the angular component decreases and 
vanishes at He (apparently linearly in He — H). The absolute value of M increases almost 
homogeneously over the sample. A surface layer of lower magnetization is preserved even in 
the presence of strong fields. Above He the increase of (My) is only due to the increase of 
|M| and will saturate in the limit H ^ 00. 

These results should be compared with those of a micromagnetic calculation for isotropic 
spheres by Aharoni and Jakubovics ||36|. They were intended as a first step in modelling 



amorphous solid ferromagnets but due to the assumption of a vanishing anisotropy constant 
they should also be of interest for liquid ferromagnets. These authors find a vortex structure 
with the axis parallel to the external field, too. In the center of the vortex the magnetization 
points into the field direction, also for if — > 0. However, for the cube this contribution to the 
net magnetization is compensated by the opposite orientation of the magnetization near the 
edges (see Fig. |(b)). This implies that for cubes (My) ~ for — > 0, whereas for spheres 
a spontaneous net mag netization (My) {H ^ 0+) = - (M||) {H ^ 0") ^ remains in the 
limits H ^ 0^. This discontinuity at if = yields an infinite zero field susceptibility. Thus 
spontaneous liquid ferromagnetism should be easier to detect by macroscopic magnetization 
measurements in spherical than in cubic samples. 

The same transition between an inhomogenous and a nearly homogeneous state arises 
when the external field is kept fixed and the temperature is varied. As shown in Fig. ^ 
the angular magnetization component decreases with increasing T and vanishes at a crit- 
ical temperature Tc{H), which is lower than Tc{H = 0) (compare Fig. |Tl|). The parallel 
component is nearly constant below Tc and decreases gradually above Tc where the sample 
is almost homogeneously magnetized. Our previous work pO| as well as Ref. [0 predict 



a kink in the curve M||(T) and the constant value Ai{T < T^) = H/^A-kD) with the de- 
magnetization factor D. In the present units, using i^ = 1/3 for the cube, this would mean 
M(T < Tc) = {3/ATT)^/^H/{pm) ~ 0.0827 for the parameters used in Fig. [l^. Both features 
have also been observed in ellipsoidal samples of a solid dipolar Ising ferromagnet [3S]. Ob- 
viously the present theory yields a higher plateau value of the magnetization and a rounding 
of the kink at the phase transition. One also observes that the angular component vanishes 
linearly at the critical point, as in Fig. 0, but in contrast to the simplified model used in 
Ref. We surmise that these differences are due to the fact that for a cube the assumption 
adopted in Refs. and |2^ of a homogeneous magnetization My in both phases is not 
fulfilled. Another possible source for this dicrepancy could be the finite size of the systems 
we examined. 



V. SUMMARY 

Based on density-functional theory we have obtained the following main results for the 
long-ranged orientational order of a dipolar liquid confined to a cube: 

1. The equilibrium magnetization configuration corresponds to a predominantly planar 
vortex structure which can be viewed as being composed of four triangular domains 
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separated by thick domain walls along the diagonal planes of the cube. Near one of 
the symmetry axes of the cube chosen spontaneously the magnetization escapes into 
the third dimension in order to avoid a topological line singularity (Figs. |l|, 0, and H). 

2. Point defects arise at the centers of the top and bottom face of the cube. The absolute 
value |M| of the magnetization is reduced inside the cores of these defects (Fig. ^ 
whose sizes scale with the system size. Therefore an appropiate description of the 
structure requires to take the spatial variation of |M| into account. 

3. Near the surfaces there is a layer of reduced orientational order (Figs. § and 

4. The domain walls (Fig. |[) are mainly of the Neel type, i.e., upon traversing the wall 
the magnetization rotates within the plane spanned by its asymptotic directions deep 



inside the adjacent domains. The size dependence of the wall thickness (Fig. p!0| ) and of 
other features of the configurations are compatible with the scaling behavior formulated 
in Eq. {Mj. 



The ferromagnetic order vanishes at a critical temperature which depends on the sys- 
tem size (Fig. |Tl|). Upon increasing the temperature the "escape" region near the 
vortex axis is gradually replaced by a column of disordered fluid (Figs. ^ and |T3p. 



In weak external fields normal to the surfaces of the cube the vortex axis is aligned 
parallel to the field direction and the overall structure is similar to the one in zero field 
(Figs. |T^ and jTS]). For stronger fields the magnetization is rotated increasingly into 



the field direction until a transition to an approximately homogeneously magnetized 
state takes place at a critical field strength beyond which the vortex structure is 
lost (Fig. 0). In an external field the magnetization components normal to the field 
direction vanish linearly upon approaching a critical temperature (Fig. |r^). 

Recently metallic liquid ferromagnets have been discovered in undercooled CoPd alloys 
For the formation of the ferromagnetic order of these materials short-ranged ex- 
change interactions play an important role. However, in addition dipolar interactions are 
also present and, as in solid ferromagnets, are essential in forming the domain structure. 
Therefore for these systems we expect strong similarities to the behavior of dipolar fiuids 
as discussed in this work. Up to now these undercooled liquid alloys have been prepared 
only as electromagnetically suspended spherical samples. One can speculate that the do- 
main structure within a sphere has a vortex axis, too, but no domain walls due to the higher 
symmetry compared to a cube. Our analysis indicates that an experimental study of these 
structures (e.g., by using magnetic neutron or X-ray tomography |^) would certainly be 
very rewarding. 
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FIGURES 



FIG. 1. Sections perpendicular to the z axis through the magnetization structure of a ferro- 
magnetic Uquid in a cubic volume (for L* = 7.2 and N = 20); (a) section plane z* = 0.18 near the 
midplane, (b) section plane z* = —3.06 near the bottom. The arrows represent projections of the 
local magnetization (at their midpoint) onto the section plane. 

FIG. 2. The same as in Fig. |l|, but sections perpendicular to the y axis; (a) y* = —1.62, (b) 
y* = —0.54. The schematic drawing in (c) demonstrates the mechanism of the "escape into the 
third dimension" which avoids the formation of a topological line defect along the z axis. 9 denotes 
the polar angle of the magnetization. (The increase of 9 to values larger than 7r/2 close to the 
sample edges (see Fig. §(b)) is not shown.) 

FIG. 3. The absolute value of the dimensionless magnetization near the point defect at the 
center of the bottom face of the cube (for L* = 7.2 and N = 24). |M| is plotted as a function of 
the distance from the z axis along lines parallel to the x or y axis for a series of fixed values of z 
(compare the case z ~ in, c.f.. Fig. ||(a)). One finds a pronounced decrease of |M| near the core 
of the point defect. 

FIG. 4. The cosine of the polar angle of the magnetization (cos 9 = M3/M) along (a) the center 
and (b) the diagonal lines within a plane close to the midplane {z^ — a/2, with the lattice constant 
a = L/N) for = 20 and different system sizes L*. r = \/ x'^ + y"^ is the distance from the center 
of the cube; r/L < 0.5 for the center lines (a) and r/L < l/\/2 = 0.707 for the diagonals (b). The 
values of cos 9 are averaged over the four equivalent sites for each value of r/L. The magnetization 
"escapes" into the positive z direction near the center {cos9 — > 1), lies in plane {cos9 = 0) near 
the sample surfaces (center lines), and turns into the negative z direction (cos 9 = — 1) near the 
edges. Except for the smallest system (L* = 4.8) the proposed scaling with the system size (see 
Eq. (|26|)) is fulfilled approximately, i.e., in the limit L ^ a two master curves evolve, one for the 
center lines and one for the diagonals. 

FIG. 5. The absolute value M of the dimensionless magnetization M for the same parameters 
as in Fig. |^ averaged over (a) the four center lines and (b) the four diagonals. With the exception 
of the smallest system M is approximately constant in the bulk and decreases near the surfaces. 
For L ^ a the behavior turns into a single master curve consistent with the proposed scaling 
(Eq. dH)). 

FIG. 6. Averaged absolute value {M)xy of the magnetization within the planes z = const for 
= 20 and several system sizes. It is constant in the inner part of the sample and decreases near 

the bottom and top surfaces. The size of the disordered region remains proportional to L for large 

L. 
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FIG. 7. Reduced free energy difference between the isotropic and the ferromagnetic phase 
fif ~ i^if / /^) ^ ^ function of the lattice constant a/a = L*/N of the discretization mesh 
for three system sizes. The finite lattice constant induces minima (maxima) near the points where 
the inverse a /a = N/L* is an integer (one half plus an integer). These oscillations die out for 
a ^ 0. 

FIG. 8. Variation of the transversal magnetization component Mt^xy = {Mi + M2)/\/2 
upon traversing the domain wall along a normal to the wall through the point 
Ro/cj = (-1.95,-1.95,-0.15) (see the inset) for L* = 7.2 and iV = 24 (compare Fig. g(a)). 
This is the component which characterizes the difference of the magnetization directions of the 
adjacent domains. Our definition of the wall thickness S can be inferred from the figure. Here 
r/L = 0.38. 

FIG. 9. Dependence of the domain wall thickness 6 on the distance r from the center of the 
plane z = —a/2 ~ (i.e., close to the midplane) for L* = 7.2 and different numbers of lattice 
points N, averaged over the four equal walls in the cube (see Fig. ^). In the medium range of r the 
thickness 6 depends only weakly on r and A^, while near the center and edges 6 is influenced by 
the escape of the orientation into the z direction and the vicinity of the surfaces, respectively. 

FIG. 10. The wall thickness at half distance between the center and the edges of the cube 
as a function of the system size L. The different points at the same L correspond to different 
lattice constants. (From top to bottom: L* = 4.8, N = 10, 12, 14, 16, 20, 18, 22, 24; L* = 7.2, 

= 14,12,18,10,22,16,20,24; L* = 9.6, N = 18,14,16,24,20,22; L* = 12, iV = 22,20.) The 
extrapolation to 1/L = suggests a finite value of 5/L in the thermodynamic limit. 

FIG. 11. The squared spatially averag ed order parameter (|M|)2 = (iV^^ |]y[(-j^)|^^ 
creases linearly with increasing temperature T*. It vanishes at a critical temperature Tc{L) which 
depends on the system size. This linear behavior holds even outside the close vicinity of Tc{L) 
where it has to be so due the inherent mean-field character of the present theory. 

FIG. 12. The polar angle of M (see Figs. §(c) and ^) along the central lines near the midplane 
(fixed z = —a/2 as in Fig. ^ as a function of the temperature for L* = 9.6 and = 16. The 
escape near the vortex axis reduces upon approaching the critical point T*(L* = 9.6) = 2.89. 

FIG. 13. The absolute value of the magnetization along the diagonal lines (for fixed z = —a/2) 
at various temperatures. The different line styles correspond to the same temperatures as in Fig. |l^. 
Near the phase transition T*{L* = 9.6) = 2.89 the orientational order is particularly reduced near 
the vortex axis and in a surface layer of increasing thickness. 
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FIG. 14. Sections orthogonal to the field direction through the magnetization structure for 
L* = 9.6 and N = 16 at z* = -0.9 (a) without and (b) with an external field H* = R^Ja^e = 1 
applied in the z direction. The scale factor which determines the lengths of the arrows is the same 
in both parts of this figure and also in Fig. 15. The transversal components of the magnetization 
are reduced by the field but the overall feature of the inhomogeneous structure is preserved. 



FIG. 15. Vertical sections perpendicular to the y axis and thus parallel to the field direction 
through the same magnetization structures as in Fig. 14 for y* = —2.1. At this distance of the 
plane from the center the z component of the magnetization in zero field is only small, while the 
external field induces a large z component everywhere. 



FIG. 16. Dependence of the averaged longitudinal magnetization component (M||) and the 
averaged angular component (M(^) on the field strength (L* = 9.6, N = 16, T* = 2.25). The 
angular component decreases for increasing and vanishes seemingly linearly at ff* ~ 1.50. Below 
H* the parallel component increases approximately linearly and crosses over towards saturation 
for a* > H*. Above H* the sample has an approximately homogeneous magnetization. Compare 
Fig. 11 in Ref. g^. 

FIG. 17. The same magnetization components as in Fig. ^ for fixed H* = 1 as a function 
of temperature (L* = 9.6, N = 16). There is a transition from a vortex structure to an almost 
homogeneous magnetization at T*{H*) ~ 2.65. Compare Fig. 12 in Ref. pO[| . 
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Fig. 12 
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